Adaptive Walks

Imports

library(notly)

Attaching package: ‘notly’

The following object is masked from ‘package:plotly’:

    ggplotly

View a fitness landscape

p <- plotFitnessLandscape()

save_dir <- file.path(output_dir, "FitnessFunctions")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
fname <- file.path(save_dir, "fitness_function.html")
htmlwidgets::saveWidget(as_widget(p), fname)

Plot different trait architectures according to different algorithms

save_dir <- file.path(output_dir, "TraitArchitecture")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

# Additive effects
methodType <- "Additive"
g <- plotTraitArchitecture(founderPop, methodType)
ggplot2::ggsave(filename = paste0("traitarchitecture_", methodType, ".pdf"),
                path=save_dir,
                device = "pdf")

# Effect of each allele on fitness
methodType <- "Fitness"
g <- plotTraitArchitecture(founderPop, methodType)
ggplot2::ggsave(filename = paste0("traitarchitecture_", methodType, ".pdf"),
                path=save_dir,
                device = "pdf")

p <- plot3dPopulationFitness(founderPop, calculateFitnessTwoTrait)

fname <- file.path(save_dir, "3dpopulationfitness.html")
htmlwidgets::saveWidget(as_widget(p), fname)

Simulate several adaptive walks with different population sizes

# Different starting population sizes
subPopSizes <- c(50,500)
selProps <- c(0.4,0.4)
n.h2 <- 0.3
n.var <- 0.1
n.qtlPerChr <- 2
n.gens <- 50

# Don't use a SNP chip for these simulations
addSnpChip <- FALSE

fig <- plot_ly() %>%
  layout(legend=list(title=list(text='Population Size')),
         scene = list(xaxis = list(title = "Trait A"),
                      yaxis = list(title = "Trait B"),
                      zaxis = list(title = "Fitness"),
                      aspectmode='cube')) #hide_colorbar()

#source("Scripts/CreateFounderPop.R")
pop <- founderPop
for (gen in 1:n.burnInGens) {
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=nInd(pop)*n.burnInSelProp, nCrosses=nInd(pop))
}

# Iterate through each population size
for (p in 1:length(subPopSizes)) {
  n.subPopSize <- subPopSizes[p]
  n.selProp <- selProps[p]
  print(n.subPopSize)
  fit_df <- data.frame(gen=1:n.gens,
                   fitness=numeric(n.gens),
                   traitValA=numeric(n.gens),
                   traitValB=numeric(n.gens))
  subPop <- selectInd(pop, nInd=n.subPopSize, use="rand")
  # Iterate through the generations, update the result dataframe, and advance
  # progeny based on the two trait fitness funcion
  for(gen in 1:n.gens) {
    meanFitness <- mean(twoTraitFitFunc(pheno(subPop)))
    fit_df$fitness[gen] <- calculateFitnessTwoTraitModified(meanP(subPop)[1], meanP(subPop)[2])
    fit_df$traitValA[gen] <- meanP(subPop)[1]
    fit_df$traitValB[gen] <- meanP(subPop)[2]
    selRat <- selectionRatio(meanFitness)
    subPop <- selectCross(subPop, trait=twoTraitFitFunc, nInd=n.subPopSize*selRat, nCrosses=n.subPopSize)
  }
  # Add a trace to represent this population's adaptive walk
  fig <- fig %>% add_trace(
    fig,
    fit_df,
    name = n.subPopSize,
    x = fit_df$traitValA,
    y = fit_df$traitValB,
    z = fit_df$fitness,
    type = 'scatter3d',
    mode = 'lines',
    opacity = 1,
    line = list(autocolorscale=FALSE, colors="PuBuGn", color=p, which=2, width = 10)
  )
}

[1] 50
[1] 500

# Create a matrix of fitness values, with a small increment along the x and y axes.
fitness_x = seq(n.initTraitVal*-1,n.initTraitVal*1, by=(n.initTraitVal/100))
fitness_y = seq(n.initTraitVal*-1,n.initTraitVal*1, by=(n.initTraitVal/100))
fitness_z = outer(fitness_x,fitness_y,calculateFitnessTwoTrait)

fig <- fig %>%
      add_trace(
        fig,
        x=fitness_x,
        y=fitness_y,
        z=fitness_z,
        type='surface',
        colorbar=list(title = "Fitness"),
        colors = viridis(n=10),
        opacity=1.0)


fig

save_dir <- file.path(output_dir, "DifferentPopulationSizes")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

fname <- file.path(save_dir, paste0("populationSizes_", subPopSizes[1], ",", subPopSizes[2], ",", subPopSizes[3], ".html"))
htmlwidgets::saveWidget(as_widget(fig), fname)

Overlay an adaptive walk over a fitness landscape, and save 3D and 2D versions

fig <- plot_ly()
fit_df <- data.frame(gen=1:n.gens,
                   fitness=numeric(n.gens),
                   traitValA=numeric(n.gens),
                   traitValB=numeric(n.gens))

pop <- founderPop

for(gen in 1:n.gens) {
  fit_df$fitness[gen] <- mean(twoTraitFitFunc(gv(pop)))
  fit_df$traitValA[gen] <- meanG(pop)[1]
  fit_df$traitValB[gen] <- meanG(pop)[2]
  meanFitness <- mean(twoTraitFitFunc(pheno(pop)))
  selRat <- selectionRatio(meanFitness)
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=n.popSize*selRat, nCrosses=n.popSize)
}

save_dir <- file.path(output_dir, "OverlayAdaptiveWalk")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)

plotType <- "CONTOUR"
pc <- overlayWalkOnLandscape(fit_df, type=plotType, calculateFitnessTwoTrait)
fname <- file.path(save_dir, paste0("adaptivewalk_", plotType, ".html"))
htmlwidgets::saveWidget(as_widget(pc), fname)

plotType <- "SURFACE"
ps <- overlayWalkOnLandscape(fit_df, type=plotType, calculateFitnessTwoTrait)
fname <- file.path(save_dir, paste0("adaptivewalk_", plotType, ".html"))
htmlwidgets::saveWidget(as_widget(ps), fname)

write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

Plot the change in allele frequency over time

n.qtlPerChr <- 2
source("Scripts/CreateFounderPop.R")
pop <- founderPop

# Results dataframe
freq.df <- data.frame(gen=1:n.gens)

# Get the effect sizes of each qtl
qtlEff.df <- getQtlEffectSizes(pop)

# Get the names of all the QTLs
qtl <- colnames(getUniqueQtl(pop))

# Create a dataframe of all zeros where the columns are the QTL ids, and the # rows is the # of generations
qtl.df <- data.frame(matrix(0, ncol=length(qtl), nrow=n.gens))
colnames(qtl.df) <- qtl

# Combine the dataframes
freq.df <- cbind(freq.df, qtl.df)

# Iterate through each generation
for(gen in 1:n.gens) {
  meanFitness <- mean(twoTraitFitFunc(pheno(pop)))
  # Get the qtl genotype data
  qtlGeno <- getUniqueQtl(pop)
  # Get the frequency of the '2' allele at each locus
  for (l in 1:length(qtl)) {
    # id is the name of the qtl (chr_site)
    id <- qtl[l]
    # A list of genotype data for each individual in the population at that locus
    locus <- qtlGeno[,id]
    # Calculate the allele frequency as the frequency of homozygous individuals (for 'allele') +
    # 1/2 * frequency of heterozygous individuals (assumes the locus is biallelic)
    freq.df[gen,id] <- (sum(locus==n.allele)/n.popSize) + ((sum(locus==1)/n.popSize)/2)
  }
  # Determine selection ration based on fitness
  selRat <- selectionRatio(meanFitness)
  #Advance individuals
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=n.popSize*selRat, nCrosses=n.popSize)
}

# Make the dataframe tidy
freq.df<- melt(freq.df, id="gen", variable.name="QTL ID", value.name="Allele Frequency")

# Add the qtl effect size data to the dataframe
freq.df <- merge(freq.df, qtlEff.df, by.x="QTL ID", by.y="row.names", all.x=TRUE)

# Create a line plot for the change in frequency of alleles over time
# Each line's opacity is a function of its effect size (higher=darker)
g <- ggplot(freq.df, aes(x=gen, y=freq, color=id, alpha=eff_size)) +
  geom_line(size=0.7)

save_dir <- file.path(output_dir, "AlleleFrequencies")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

ggplot2::ggsave(filename = paste0("allelefrequencies.pdf"),
                path=save_dir,
                device = "pdf",
                width=10,
                height=7)

Simulate an adaptive walk, and plot the following curve: #segregating QTL/#segregating loci

Original question: out of segregating alleles left in the population, how many of them would increase fitness? Is this even a valid question?

TODO: - How does this change with size of allele? - How does size of ‘fitness delta’ change over generations? - Simualate mutations for a single individual to figure out P(allelic substitution is favorable), should reflect FKO - P(mutation gets fixed) - show that this matches Kimura - Compare adaptive walks for DH/inbred population where there are new mutations VS population w/standing genetic variation

n.gens <- 200

# Reset variables
source("Scripts/CreateFounderPop.R")
source("Scripts/GlobalParameters.R")

pop <- founderPop
df <- data.frame(gen=c(),
                 fitness=c(),
                 traitVal1=c(),
                 traitVal2=c(),
                 percFitInc=c(),
                 segAlleles=c(),
                 segQtl=c())

# Iterate through the generations
for (gen in 1:n.gens) {
  print(gen)
  # Counter variables
  numHetLoci = 0
  numHetQtl = 0
  numInc = 0

  # Get all of the loci from the population
  segSiteGeno <- pullSegSiteGeno(pop)
  # Get all of the QTL from the population (a sampling of the segSiteGeno)
  qtlGeno <- getUniqueQtl(pop)
  qtlLoci <- colnames(qtlGeno)
  loci <- colnames(segSiteGeno)
  nLoci <- ncol(segSiteGeno)
  # Iterate through all of the loci
  # TODO: just calculate numHetQtl / numHetLoci
  for (l in 1:ncol(segSiteGeno)) {
    id <- loci[l]
    locus = segSiteGeno[,l]
    # Check if the locus is segregating
    if (hetLocus(locus)) {
      # Increment the counter
      numHetLoci <- numHetLoci + 1
      # Determine the effect size - this is where the bug is
      #e <- getEffectSize(locus, id, pop, "Fitness")
      #if (e > 0) {
      #  numInc <- numInc +1
      #}
      # If the locus is a QTL, then it has an effect size
      if (id %in% qtlLoci) {
        numHetQtl <- numHetQtl + 1
      }
    }
  }

  # Calculate the ratio of segregating QTL to segregating alleles
  if (numHetLoci == 0) {
    perc <- 0
  } else {
    perc <- (numHetQtl / numHetLoci)
  }
  df <- rbind(df, data.frame(gen=gen,
                             fitness=mean(twoTraitFitFunc(gv(pop))),
                             traitVal1=meanG(pop)[1],
                             traitVal2=meanG(pop)[2],
                             percFitInc=perc,
                             segAlleles=numHetLoci,
                             segQtl=numHetQtl))
  # If the population is within n.margin, terminate the simulation
  if (mean(twoTraitFitFunc(gv(pop))) >= n.margin) {
    break
  }
  meanFitness <- mean(twoTraitFitFunc(pheno(pop)))
  selRat <- selectionRatio(meanFitness)
  # Advance the top progeny according to the fitness function
  pop <- selectCross(pop=pop, trait=twoTraitFitFunc, nInd=nInd(pop)*selRat, nCrosses=nInd(pop))
}

save_dir <- file.path(output_dir, "BeneficialAlleles")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

fname <- file.path(save_dir, "beneficialalleles.pdf")
pdf(fname)

# Plot the ratio of segregating QTL to segregating alleles, along with fitness
par(mar = c(5, 5, 3, 5) + 0.3)
plot(df$gen, df$fitness, type="l", lwd = "3", col=2, xlab = "Generation", ylab = "Fitness")
par(new = TRUE) 
plot(df$gen, df$percFitInc, type="l", lwd = "3", col = 3, axes = FALSE, xlab = "", ylab = "") 
axis(side = 4, at = pretty(range(df$percFitInc)))
mtext("% of Segregating alleles that are QTL", side = 4, line = 3)
par(xpd=TRUE)
legend("right",
  c("Fitness", "% Alleles"),
       lty = 1,
       lwd = 3,
       col = 2:3)
dev.off()

This block runs n.nPops * n.sims Monte Carlo simulations of different populations to determine the change in the average effect size of alleles that are fixed along the adaptive walk, saving the order in which the alleles are fixed n.nPops populations are created, and each one undergoes n.sims unique adaptive walks

This block will simulate one base population, from which two sub-populations are selected, and undergo purifying selection independently.

source("Scripts/GlobalParameters.R")
source("Scripts/CreateFounderPop.R")

fit_df <- data.frame(gen=1:n.burnInGens,
                 fitness=numeric(n.burnInGens),
                 traitValA=numeric(n.burnInGens),
                 traitValB=numeric(n.burnInGens))
pop <- founderPop

# Burn-in generations
for (gen in 1:n.burnInGens) {
  fit_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(pop)))
  fit_df$traitValA[gen] <- meanP(pop)[1]
  fit_df$traitValB[gen] <- meanP(pop)[2]
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=nInd(pop)*n.burnInSelProp, nCrosses=nInd(pop))
}

# Create a random vector of size n.pops, with a random order of sub-population ids
randVec <- sample(rep(c(1:n.nPops), times=n.popSize/n.nPops))

# Select all of the "1" indexed individuals
popA <- selectInd(pop, trait=selectSubPop, selectTop=TRUE, nInd=n.subPopSize, idx=1, randVec=randVec)
# Select all of the "2" indexed individuals
popB <- selectInd(pop, trait=selectSubPop, selectTop=TRUE, nInd=n.subPopSize, idx=2, randVec=randVec)

# Create dataframes for each subpopulation, initializing with current values
popA_df <- data.frame(gen=c(1),
                 fitness=c(mean(twoTraitFitFunc(pheno(popA)))),
                 traitValA=c(meanP(popA)[1]),
                 traitValB=c(meanP(popA)[2]))
popB_df <- data.frame(gen=c(1),
                 fitness=c(mean(twoTraitFitFunc(pheno(popB)))),
                 traitValA=c(meanP(popB)[1]),
                 traitValB=c(meanP(popB)[2]))

# Iterate through the generations
for (gen in 1:n.gens) {
  # If popA is within the margin of the fitness optimum, don't progress it any further
  if (mean(twoTraitFitFunc(pheno(popA))) < n.margin) {
    # Advance the population
    meanFitness <- mean(twoTraitFitFunc(pheno(popA)))
    # Get a selection ratio based on fitness
    selRat <- selectionRatio(meanFitness)
    popA <- selectCross(popA, trait=twoTraitFitFunc, nInd=nInd(popA)*selRat, nCrosses=nInd(popA))
    # Update the dataframe with new values
    popA_df <- rbind(popA_df, data.frame(gen=gen,
                                       fitness=meanFitness,
                                       traitValA=meanP(popA)[1],
                                       traitValB=meanP(popA)[2]))
    
  }
  # If popB is within the margin of the fitness optimum, don't progress it any further
  if (mean(twoTraitFitFunc(pheno(popB))) < n.margin) {
    # If popA is within the margin of the fitness optimum, don't progress it any further
    meanFitness <- mean(twoTraitFitFunc(pheno(popB)))
    # Get a selection ratio based on fitnes
    selRat <- selectionRatio(meanFitness)
    popB <- selectCross(popB, trait=twoTraitFitFunc, nInd=nInd(popB)*selRat, nCrosses=nInd(popB))
    # Update the dataframe with new values
    popB_df <- rbind(popB_df, data.frame(gen=gen,
                                       fitness=meanFitness,
                                       traitValA=meanP(popB)[1],
                                       traitValB=meanP(popB)[2]))
  }
}
# Update rownames
rownames(popA_df) <- 1:nrow(popA_df)
rownames(popB_df) <- 1:nrow(popB_df)

# Plot the adaptive walks
fig <- plot_ly()
fig <- add_trace(
  fig,
  fit_df,
  name = "Burn In",
  x = fit_df$traitValA,
  y = fit_df$traitValB,
  z = fit_df$fitness,
  type = 'scatter3d',
  mode = 'lines',
  opacity = 1,
  color = 'yellow',
  line = list(width = 5)
)

fig <- add_trace(
    fig,
    popA_df,
    name = "Pop A",
    x = popA_df$traitValA,
    y = popA_df$traitValB,
    z = popA_df$fitness,
    type = 'scatter3d',
    mode = 'lines',
    opacity = 1,
    color = 'red',
    line = list(width = 5)
  )

fig <- add_trace(
    fig,
    popB_df,
    name = "Pop B",
    x = popB_df$traitValA,
    y = popB_df$traitValB,
    z = popB_df$fitness,
    type = 'scatter3d',
    mode = 'lines',
    opacity = 1,
    color = 'blue',
    line = list(width = 5)
  )


p <- fig %>%
  layout(legend=list(title=list(text='Population')),
         showlegend=FALSE,
         scene = list(xaxis = list(title = "Trait A"),
                      yaxis = list(title = "Trait B"),
                      zaxis = list(title = "Fitness"),
                      aspectmode='cube')) %>% hide_colorbar()

save_dir <- file.path(output_dir, "DivergingPopulations")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
fname <- file.path(save_dir, "adaptivewalks.html")
htmlwidgets::saveWidget(as_widget(p), fname)

fname <- file.path(save_dir, "2PopulationFitness.html")
p <- plot3dPopulationFitnessTwoPops(popA, popB)
htmlwidgets::saveWidget(as_widget(p), fname)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

fname <- file.path(save_dir, "traitarchitecture.pdf")
pdf(fname)

p1 <- plotTraitArchitecture(popA, "Fitness", "popA")
p2 <- plotTraitArchitecture(popB, "Fitness", "popB")
p3 <- plotTraitArchitecture(popA, "Additive", "popA")
p4 <- plotTraitArchitecture(popB, "Additive", "popB")

(p1|p2)/(p3|p4)
dev.off()

# Create a density plot of trait 1
trait1.df <- as.data.frame(cbind(pheno(popA)[,1], pheno(popB)[,1]))
colnames(trait1.df) <- c("popA", "popB")
trait1.df <- trait1.df %>%
  pivot_longer(c("popA", "popB"), names_to="pop", values_to="pheno")
t1 <- ggplot(trait1.df, aes(pheno, fill=pop, color=pop)) +
  geom_density(alpha=0.1) +
  labs(title="Trait 1")

# Create a density plot of trait 2
trait2.df <- as.data.frame(cbind(pheno(popA)[,2], pheno(popB)[,2]))
colnames(trait2.df) <- c("popA", "popB")
trait2.df <- trait2.df %>%
  pivot_longer(c("popA", "popB"), names_to="pop", values_to="pheno")
t2 <- ggplot(trait2.df, aes(pheno, fill=pop, color=pop)) +
  geom_density(alpha=0.1) +
  labs(title="Trait 2")


(t1|t2)
ggplot2::ggsave(filename = "trait_distributions.pdf",
                path=save_dir,
                device = "pdf")

source("Scripts/GlobalParameters.R")
n.qtlPerChr <- 2
n.h2 <- 0.3
n.qtlPerChr <- 2
n.gens <- 100
#source("Scripts/CreateFounderPop.R")
n.selProp <- 0.5
n.subPopSize <- 50

fit_df <- data.frame(gen=1:n.burnInGens,
                 fitness=numeric(n.burnInGens),
                 traitValA=numeric(n.burnInGens),
                 traitValB=numeric(n.burnInGens))
pop <- founderPop

# Burn-in generations
for (gen in 1:n.burnInGens) {
  fit_df$fitness[gen] <- mean(twoTraitFitFunc(pheno(pop)))
  fit_df$traitValA[gen] <- meanP(pop)[1]
  fit_df$traitValB[gen] <- meanP(pop)[2]
  pop <- selectCross(pop, trait=twoTraitFitFunc, nInd=nInd(pop)*n.burnInSelProp, nCrosses=nInd(pop))
}

# Create a random vector of size n.pops, with a random order of sub-population ids
randVec <- sample(rep(c(1:n.nPops), times=n.popSize/n.nPops))

# Select all of the "1" indexed individuals
popA <- selectInd(pop, trait=selectSubPop, selectTop=TRUE, nInd=n.subPopSize, idx=1, randVec=randVec)
# Select all of the "2" indexed individuals
popB <- selectInd(pop, trait=selectSubPop, selectTop=TRUE, nInd=n.subPopSize, idx=2, randVec=randVec)

# Create dataframes for each subpopulation, initializing with current values
popA_df <- data.frame(gen=c(1),
                 fitness=c(mean(twoTraitFitFunc(pheno(popA)))),
                 traitValA=c(meanP(popA)[1]),
                 traitValB=c(meanP(popA)[2]))
popB_df <- data.frame(gen=c(1),
                 fitness=c(mean(twoTraitFitFunc(pheno(popB)))),
                 traitValA=c(meanP(popB)[1]),
                 traitValB=c(meanP(popB)[2]))

# Iterate through the generations
for (gen in 1:n.gens) {
  # If popA is within the margin of the fitness optimum, don't progress it any further
  if (mean(twoTraitFitFunc(pheno(popA))) < n.margin) {
    # Advance the population
    meanFitness <- mean(twoTraitFitFunc(pheno(popA)))
    # Get a selection ratio based on fitness
    selRat <- selectionRatio(meanFitness)
    popA <- selectCross(popA, trait=twoTraitFitFunc, nInd=nInd(popA)*selRat, nCrosses=nInd(popA))
    # Update the dataframe with new values
    popA_df <- rbind(popA_df, data.frame(gen=gen,
                                       fitness=meanFitness,
                                       traitValA=meanP(popA)[1],
                                       traitValB=meanP(popA)[2]))
    
  }
  # If popB is within the margin of the fitness optimum, don't progress it any further
  if (mean(twoTraitFitFunc(pheno(popB))) < n.margin) {
    # If popA is within the margin of the fitness optimum, don't progress it any further
    meanFitness <- mean(twoTraitFitFunc(pheno(popB)))
    # Get a selection ratio based on fitnes
    selRat <- selectionRatio(meanFitness)
    popB <- selectCross(popB, trait=twoTraitFitFunc, nInd=nInd(popB)*selRat, nCrosses=nInd(popB))
    # Update the dataframe with new values
    popB_df <- rbind(popB_df, data.frame(gen=gen,
                                       fitness=meanFitness,
                                       traitValA=meanP(popB)[1],
                                       traitValB=meanP(popB)[2]))
  }
}
# Update rownames
rownames(popA_df) <- 1:nrow(popA_df)
rownames(popB_df) <- 1:nrow(popB_df)

# Plot the adaptive walks
# Create a matrix of fitness values, with a small increment along the x and y axes.
fitness_x = seq(n.initTraitVal*-1,n.initTraitVal*1, by=(n.initTraitVal/100))
fitness_y = seq(n.initTraitVal*-1,n.initTraitVal*1, by=(n.initTraitVal/100))
fitness_z = outer(fitness_x,fitness_y,calculateFitnessTwoTrait)

fig <- plot_ly() %>%
      layout(xaxis = list(title = "Trait A", constrain = "domain"),
             yaxis = list(title = "Trait B", scaleanchor="x")) %>%
      add_trace(
        fig,
        x=fitness_x,
        y=fitness_y,
        z=fitness_z,
        type='contour',
        colors = viridis(n=10),
        colorbar=list(title = "Fitness"),
        line = list(color = 'black', width = 1),
        opacity=1) %>%
      add_trace(
        fig,
        fit_df,
        name = "Burn-in",
        x = fit_df$traitValA,
        y = fit_df$traitValB,
        type='scatter',
        mode = 'lines',
        line = list(color = 'orange', width = 4, dash = 'solid'),
        opacity = 1) %>%
      add_trace(
        fig,
        popA_df,
        name = "Subpop 1",
        x = popA_df$traitValA,
        y = popA_df$traitValB,
        type='scatter',
        mode = 'lines',
        line = list(color = 'red', width = 4, dash = 'solid'),
        opacity = 1) %>%
      add_trace(
        fig,
        popB_df,
        name = "Subpop 2",
        x = popB_df$traitValA,
        y = popB_df$traitValB,
        type='scatter',
        mode = 'lines',
        line = list(color = 'blue', width = 4, dash = 'solid'),
        opacity = 1)

save_dir <- file.path(output_dir, "DivergingPopulations")
if (!dir.exists(save_dir)) dir.create(save_dir)
save_dir <- file.path(save_dir, format(Sys.time(), "%F_%H_%M_%S"))
if (!dir.exists(save_dir)) dir.create(save_dir)
fname <- file.path(save_dir, "adaptivewalks.html")
htmlwidgets::saveWidget(as_widget(fig), fname)

fname <- file.path(save_dir, "2PopulationFitness.html")
p <- plot3dPopulationFitnessTwoPops(popA, popB)
htmlwidgets::saveWidget(as_widget(p), fname)
write.table(getParams(), file.path(save_dir, "params.txt"), col.names=FALSE, quote=FALSE, sep=":\t")

fname <- file.path(save_dir, "traitarchitecture.pdf")
pdf(fname)

p1 <- plotTraitArchitecture(popA, "Fitness", "popA")
p2 <- plotTraitArchitecture(popB, "Fitness", "popB")
p3 <- plotTraitArchitecture(popA, "Additive", "popA")
p4 <- plotTraitArchitecture(popB, "Additive", "popB")

(p1|p2)/(p3|p4)
dev.off()

null device 
          1 

# Create a density plot of trait 1
trait1.df <- as.data.frame(cbind(pheno(popA)[,1], pheno(popB)[,1]))
colnames(trait1.df) <- c("popA", "popB")
trait1.df <- trait1.df %>%
  pivot_longer(c("popA", "popB"), names_to="pop", values_to="pheno")
t1 <- ggplot(trait1.df, aes(pheno, fill=pop, color=pop)) +
  geom_density(alpha=0.1) +
  labs(title="Trait 1")

# Create a density plot of trait 2
trait2.df <- as.data.frame(cbind(pheno(popA)[,2], pheno(popB)[,2]))
colnames(trait2.df) <- c("popA", "popB")
trait2.df <- trait2.df %>%
  pivot_longer(c("popA", "popB"), names_to="pop", values_to="pheno")
t2 <- ggplot(trait2.df, aes(pheno, fill=pop, color=pop)) +
  geom_density(alpha=0.1) +
  labs(title="Trait 2")


(t1|t2)

ggplot2::ggsave(filename = "trait_distributions.pdf",
                path=save_dir,
                device = "pdf")

Saving 7.29 x 4.51 in image